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Configuration-interaction Monte Carlo method and its application to the trapped 
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We develop a quantum Monte Carlo method to estimate the ground-state energy of a fermionic 
many-particle system in the configuration-interaction shell model approach. The fermionic sign 
problem is circumvented by using a guiding wave function in Fock space. The method provides an 
upper bound on the ground-state energy whose tightness depends on the choice of the guiding wave 
function. We argue that the antisymmetric geminal product class of wave functions is a good choice 
for guiding wave functions. We demonstrate our method for the trapped two-species fermionic cold 
atom system in the unitary regime of infinite scattering length using the particle-number projected 
Hartree-Fock-Bogoliubov wave function as the guiding wave function. We estimate the ground- 
state energy and energy-staggering pairing gap as a function of the number of particles. Our results 
compare favorably with exact numerical diagonalization results and with previous coordinate-space 
Monte Carlo calculations. 
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The configuration-interaction (CI) shell model ap- 
proach is widely used in nuclear, atomic and molecular 
physics. It accounts for both shell effects and correlations 
in finite-size many-particle systems (see, e.g., Refs. [l|, @)- 
The CI many-particle model space for an A-fermion sys- 
tem is a truncated space spanned by A-particle Slater de- 
terminants constructed from a given finite single-particle 
basis. An effective CI Hamiltonian is defined in this trun- 
cated space. 

When the many-particle model space is sufficiently 
small, the CI Hamiltonian can be diagonalized by con- 
ventional methods. However, the combinatorial increase 
of the dimensionality of the many-particle space versus 
the size of the single-particle basis and/or the number of 
particles prohibits direct diagonalization in large spaces. 
This difficulty can be overcome in part by using quan- 
tum Monte Carlo methods for which the computational 
effort scales much more gently with the size of the single- 
particle model space. An example is the auxiliary-field 
Monte Carlo (AFMC) method. The AFMC approach 
has been applied within the CI framework to nuclei [3|-|6[ 
(where the method is known as the shell model Monte 
Carlo method), and more recently to the cold atomic 
condensate in a harmonic trap [7|. 

Other quantum Monte Carlo methods include 
coordinate-space ('r-space') Monte Carlo methods, e.g., 
diffusion Monte Carlo and Green's function Monte Carlo 
(GFMC) Q. These methods filter out the ground state 
with the help of an appropriately defined projection op- 
erator or a propagator. A GFMC method that works on 
a discrete lattice in coordinate space (known as lattice 
GFMC) was also developed jUTJ. 

Here we introduce a novel configuration-interaction 
Monte Carlo (CIMC) method that provides an upper 
bound of the ground-state energy of a fermionic system in 
the CI framework. Our method builds on techniques de- 
veloped in the lattice GFMC method. A discrete configu- 
ration space is defined by the occupation numbers of the 



single-particle states and an initial configuration-space 
wave function is propagated in imaginary time. The sign 
problem is circumvented by introducing a propagator 
that depends on a guiding wave function and keeps the 
wave function positive semidefinite. The method yields 
an upper bound on the ground-state energy whose accu- 
racy (as an estimate of the ground-state energy) depends 
on the choice of the guiding wave function. We argue that 
the use of antisymmetric geminal product (AGP) wave 
functions [12j | as guiding wave functions offer a good com- 
promise between accuracy and computational efficiency. 

We demonstrate the CIMC method by its application 
to the trapped two-species fermionic cold atom system 
in the unitary limit of infinite scattering length. We use 
the particle-number projected Hartree-Fock-Bogoliubov 
(HFB) wave function (a member of the AGP class) as the 
parameter-free guiding wave function. Cold atomic gases 
have recently attracted much interest both experimen- 
tally and theoretically [13, [lj] . The interaction strength 
in these systems can be controlled by tuning the scatter- 
ing length near a Fano-Feshbach resonance. Since these 
systems depend on a few number of parameters, they are 
useful for testing many-body methods of strongly inter- 
acting systems. The unitary limit is particularly chal- 
lenging since, in the absence of a small parameter, it is 
not amenable to perturbative treatments. 

Our method consists of two main components: (i) pro- 
jecting on the ground-state wave function through a ran- 
dom walk in configuration space, and (ii) circumventing 
the fermionic sign problem with the help of a guiding 
wave function. 

We assume a general CI Hamiltonian which includes 
only two-body interactions 
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where a\ creates a particle in the single-particle state la- 
beled by a. The set S of single-particle states is assumed 



to be finite of size Af s . 

We define an operator V by 



V = 1- At(H - E T ) , 



(2) 



where Ej~ is an energy shift (to be discussed later). This 
operator can be used to propagate in imaginary time a 
wave function in the many-particle space from r to t+At 
by 



|*t+At> - V |tf T ) 



(3) 



The ground-state wave function |\&g S ) can be projected 
out by the repeated application of V on an initial wave 
function ^o that has a non-zero overlap with ^ gs , i.e., 
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provided the eigenvalues of V are between —1 and 1. 
The latter condition implies that the time step At sat- 
isfies At < 2/(E max — Et), where -Einax is the maximal 
eigenvalue of H. Equation ([JJ is exact and there is no 
error that depends on the size of the time step. 

The method works in the TV-particle Hilbert space that 
is spanned by the set of all A-particle Slater determinants 
constructed from the single-particle orbitals a G S. We 
will denote these Slater determinants or 'configurations' 
by |n), where n = {n a } and n a = 0, 1 is the occupation 
number of orbital a. The one-step propagation ([3]) can 
be written in the configuration representation as 



* r+ Ar(m) 
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where ^r(n) = (n|\l/ T ) is the wave function representa- 
tion in configuration space. We rewrite the configuration- 
space matrix elements of V in the form 



where 



and 



(m|7>|n) =5(n)p(m,n) 



g(n) = ]T(m|P|n) 



p(m,n) = 



#( n ) 
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We first consider the case when the matrix elements 
of V are all non-negative, i.e., (m|'P|n) > for all m 
and n. Then < p(m, n) < 1 with ^] m p(m,n) = 1 
and <?(n) > 0, and we can interpret p(m, n) for fixed 
n as a probability and p(n) as a weight. This enables 
us to carry out the propagation of $> T in §5§ stochasti- 
cally as follows. Assuming at a given time t, the wave 
function ^ r is non- negative in configuration space, i.e., 
^r( n ) > for all n, the wave function \I> T can be repre- 
sented as an ensemble of configurations n. According to 
([5]) and the non-negativity of the matrix elements of V, 



the wave function remains non- negative at r + At, i.e., 
^r+Art 111 ) > for all m. For each configuration n, a 
new configuration m is chosen with probability p(m, n) 
and replicated with weight g(n). The resulting ensemble 
of configurations {m} samples the wave function <f r+ A r 
at the next imaginary-time step r + At. We note that 
in CIMC, the diffusion in configuration space is deter- 
mined by the interaction part of the Hamiltonian, while 
in coordinate-space Monte Carlo methods it is the kinetic 
part which governs the diffusion. 

After a sufficiently large number of time steps, the con- 
tribution from excited states to ^> T becomes negligible. 
Ensembles generated at subsequent time steps are consid- 
ered as representatives of |^ gs ) and a decorrelated subset 
of these ensembles is used to calculate the observables. 

As mentioned above, the propagation with V filters 
out the ground state when the time step satisfies At < 
2/(-E max — Et)- This upper bound for At becomes 
smaller with increasing N and/or Af s since E max gets 
larger. Consequently, the number of time steps required 
for decorrelation and for projecting the ground state 
increases. This makes the simple algorithm described 
above inefficient for large Af s and/or N. 

We overcome this latter difficulty by using an algo- 
rithm proposed in Ref. [jja. We start by choosing a fi- 
nite imaginary time step St (which is different from At) . 
Assuming the configuration n is a member of the en- 
semble at time t, we describe the choice of the corre- 
sponding configuration and its replication weight at time 
t + St. We define a propagation time St p that is ini- 
tially set to St and a weight g that is initially set to 
1. We then sample a time STd for an off-diagonal move 
from the (unnormalized) probability distribution e _7rd d , 
where 7r d = J2 m ^n( m \ H _ E T \n). If 5T d > St p , then n 
remains unchanged in the next time step and is repli- 
cated with the weight g = e-^Em^lfl-^rH. Other- 
wise, a new configuration n' is chosen with probability 
(n' \H — Ex\n) / nd (n' ^ n), and the weight factor is mul- 
tiplied by e~' 5T ' i £-< m l K- - E!T l n >. This process is repeated 
with the replacements n — > n' and St p — » 5t p — St^ until 
we have Srd > St p . This algorithm generates the prob- 
ability distribution p(m, n) and weight <?(n) that corre- 
spond to the propagator e -Sr(H-E T ) _ i} mAT ^ -pSr/Ar 

without any time step errors due to the finiteness of 
St E3. 

The above algorithm breaks down if p(m, n) < 0, i.e., 
if (m|if|n) > for any pair of configurations m ^ n. 
For a general CI Hamiltonian, this might be the case and 
is the manifestation of the Monte Carlo sign problem for 
our method 1 

In continuum r-space Monte Carlo methods the sign 
problem is circumvented with the help of a fixed-node 



1 An important exception is the pairing Hamiltonian, for which a 
diffusion Monte Carlo algorithm free of a sign problem can be 
formulated within the CI framework using the occupations of 
time-reversed pairs Il6t ll7(l . 



approximation, which can be used to obtain an upper 
bound on the ground-state energy. However, in the CI 
approach the Hilbert space is labeled by discrete quan- 
tum numbers. It was shown in lattice GFMC that a 
fixed-node-like approximation can be introduced to ob- 
tain an upper bound on the ground-state energy also for 
discrete Hilbert spaces. In the following we show how 
this can be formulated for a CI Hamiltonian. 

We choose a guiding wave function $ G and define for 
any configurations n and m with $G(n) 7^ the quantity 
s(m, n) = $G(m)(m|i?|n)/$G(n). We then introduce a 
family of Hamiltonians T-i 1 defined over configurations n 
with $G( n ) 7^ such that the off-diagonal matrix ele- 
ments are given by [9J-lll| 

H*» = { -'KK * ( " h li S ° - <W 

while the diagonal matrix elements are 
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(n|#|n) + (l+ 7 ) Yj *(m,n). (10) 



s(m,n)>0 



We note that T-L 1= -i = H. We define a 7-dependent 
propagator V-, for configurations n with $G(n) ^ by 

(m|7>» = l-Ar$ G (m)(m|H 7 -S T |n)/$ G (n) . (11) 

For 7 > the propagator V 1 satisfies (m^-yln) > 
and is therefore free of the sign problem. The 
stochastic projection algorithm outlined above can then 
be generalized with the replacement of (m|if|n) by 
<&G( m )(m|% 7 |n)/<I>G(ii). This stochastic projection fil- 
ters out the wave function $ G ( n ) v I' j( n ) , where ^ 7 (n) is 
the ground-state wave function of % 7 . 

The ground-state energy £ 7 of Hj (in the non-singular 
space <&G( n ) 7^ 0) is an upper bound for the ground-state 
energy Eq of the original Hamiltonian H for 7 > 0. This 
result is derived in Refs.llOlandllllfor the case when $g is 
nonzero for all configurations in the Hilbert space. In the 
following we derive this result for the case where $G( n ) 
may vanish for certain configurations (this is typically 
the situation for CI guiding wave functions) . In this case 
the propagator "P 7 is non-singular only in the subspace of 
configurations n spanned by $G( n ) 7^ 0. The Hamilto- 
nian 'H 7 is also non-singular in this subspace. The energy 
E 1 is an upper bound for the ground-state energy of H 
when restricted to this non-singular subspace [id lllj . 
This latter energy is an upper bound on Eq (which is 
the ground-state energy of H in the full Hilbert space). 
Thus, £ 7 is an upper bound on Eq. 

It is seen from Eq. (TTTJ) that as long as the initial 
ensemble includes only those configurations for which 
<S>g(h) 7^ 0, the Monte Carlo projection does not reach 
configurations m for which $G( m ) = 0. Thus, the Monte 
Carlo projection correctly finds the ground-state energy 
in the space where V 1 is non-singular and thus provides 
an upper bound on Eq. 



One can verify using Eqs. (j9|) and (jlQI) that 
($g|H 7 |$g) = ($g\H\$g)- Since, £ 7 is the ground- 
state energy of H 7 , £ 1 < ($g|H 7 |$g) = {$g\H\$g), 
i.e., the upper bounds on Eq given by £ 7 are tighter than 
the variational upper bound with $g- 

The energies £ 1 are estimated using the 'mixed' esti- 
mate for the Hamiltonian H 7 

££ L (n)<Mn)* 7 (n) 

r =-^ «7TrE fL ( n *)' ( 12 ) 



^d> G (n)* 7 (n) 



where Ne is the size of the ensemble {n^} represent- 
ing $ G (n)* 7 (n) and £ L (n) = ($ G |^ 7 |n)/$ G (n) = 
($ G |7J|n)/$ G (n) is the so-called local estimate of the 
energy. 

The energy shift Et is used to control the size of the 
configuration population. It is adjusted infrequently dur- 
ing the evolution to keep the population of configurations 
roughly constant. It provides an independent estimate 
(the 'growth' estimate) of the ground-state energy of H 7 . 
However, its statistical error is typically much larger than 
that of the 'mixed' estimate described above. 

The finite time step St is kept fixed throughout the 
entire calculation. In principle, its value is arbitrary since 
the method is free from time step errors. However, if St 
is too large, then the configuration population undergoes 
large fluctuations between consecutive time steps. On 
the other hand, if St is too small, then the number of 
time steps required to decorrelate the energy becomes 
very large. As a compromise we choose an intermediate 
value for St. 

A linear extrapolation of £ 1 from any two values of 
7>0to7 = — 1 also provides an upper bound on Eq 
that is tighter than the individual £ 7 's [18(]. We found 
that the best compromise between the tightness of the 
upper bound and the size of the statistical error in the 
extrapolation is obtained when the two values of 7 are 
chosen to be and 1. This choice gives -Ecimc = 2£ 7=0 — 
£ 1= i as our best upper bound for the ground-state energy 

Eq. 

The tightness of the energy upper bound is determined 
to a large extent by the quality of the guiding wave func- 
tion. If the guiding wave function is the exact ground- 
state wave function (i.e., <I> G = \& gs ), then £ 1 — Eq for all 
7. Furthermore, for $ G close to VPgs the difference £ 1 — Eq 
is second order in <I> G — VPgs [10(. Thus, the accuracy of 
the CIMC method can be improved by systematically 
improving the quality of the guiding wave function. 

Ideally, $ G should encode all our apriori knowledge 
about the ground-state wave function. However, for the 
CIMC method to be efficient, we should be able to calcu- 
late <& G quickly and accurately. As we discussed above, 
the configuration space in which the stochastic projection 
is carried out is determined by the condition $ G 7^ 0. It 
is preferable for this space to be sufficiently large so as 
to include the dominant correlations. 



In r-space Monte Carlo methods, optimized Slater- 
Jastrow or BCS-Jastrow wave functions are examples of 
highly accurate yet efficient guiding wave functions. For 
a CI Hamiltonian, a good choice is provided by the AGP 
class of wave functions [12j ■ The most general AGP wave 
function for an even ./V-particle system is given by 



|$AGP>= (X^ ofea « a 6 



N/2 



|0> 



(13) 
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where |0) is the particle vacuum. The coefficients <p a j, 
are the 'geminals' which encode information about the 
correlations in the system. The AGP wave functions in- 



corporate important two-body correlations yet they are 
described by a single pfaffian, i.e., $AGp( n ) is the pfaf- 
fian of an N x N matrix [19| . This later property ensures 
their efficient numerical evaluation. More recently, AGP 
wave functions have also been used extensively in r-space 
quantum Monte Carlo calculations (see, e.g., in Ref.l20h. 

To demonstrate the CIMC method, we consider the 
two-species (labeled by spin up and spin down) fermionic 
cold atom system, in which atoms with opposite spins in- 
teract via a short-range interaction, modeled by a contact 
interaction S(r — r'). The CI many-body Hamiltonian of 
this system in an isotropic harmonic trap is given by 



J 



H = 



J2 £ i( a n 
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where w, m and £i are, respectively, the frequency of the 
trap, the particle mass and the single-particle energies in 
an isotropic harmonic trap. The label i denotes a single- 
particle state in orbital space and f, 4- denote spin-up and 
spin-down particles. 

The dimensionless coupling strength g for finite Af 3 
is determined by the condition that the two-particle 
ground-state energy in the laboratory frame reproduces 
the exact energy 7] . In the unitary limit of infinite scat- 
tering leng th the exact two-particle ground-state energy 
is 2huj [21| . We note that the renormalization of the con- 
tact interaction in a finite model space is a non-trivial 
problem and there exist more rigorous treatments of the 
effective interaction I23-I24I. 



We choose the single-particle model space <S to include 
all single-particle orbitals within A/" ma x harmonic oscilla- 
tor shells, i.e, with energy £, < (A/" max + |)fiw. A model 
space with Af max — 9 has Af s — 440 single-particle states. 
The respective many-particle space for N — 20 particles 
has dimension of ~ 10 35 . 



For this system it is reasonable to assume that the 
dominant correlations are in the pairing and particle-hole 
channels. A wave function that includes these correla- 
tions is described by the HFB approximation. The HFB 
wave function does not conserve particle number, thus 
it is necessary to project it on a fixed number of parti- 
cles. We only consider here even-iV spin-balanced sys- 
tems and odd- A systems with one unbalanced spin up 
particle. The particle-number projected Hartree-Fock- 
Bogoliubov (PHFB) wave function (projection after vari- 
ation) belongs to the AGP class of wave functions. For 



the even-iV (N = 2n) system it is given by 



where the matrix D describes the transformation to the 
canonical basis, and (uk,Vk) are the coefficients of the 
Bogoliubov transformation in the canonical basis [2a ]. 

For odd- TV (N — 2n + 1) system, the wave function 
belongs to the generalized AGP class of wave functions 
and is given by 



I^phfb) 







(16) 
where b is the 'blocked' spin-up orbital in the canon- 
ical basis. The configuration-space representations of 
$p HFB (n) and $p HFB (n) are determinants ofnxn and 
(n + 1) x (n + 1) matrices, respectively. 

The particle-number projected Bardeen-Cooper- 
Schrieffer (PBCS) and Hartree-Fock (HF) wave func- 
tions have similar forms. In PBCS the D matrix is 
the identity, while in HF Vk/uk = 9(kp — \k\), with 
kp being the highest occupied orbital in the canonical 
basis. Generally, the quantities Dy, Uk and Vk have 
to be recalculated self-consistently depending on which 
wave function we use, e.g., the Uk,Vk in PHFB are not 
the same as those in PBCS. 

The PBCS wave function includes correlations in the 
pairing channel only, while the HF wave function includes 
correlations in the particle-hole channel only. We expect 
the PHFB wave function, which includes correlations in 
both channels to be superior to either of those. This 
expectation is confirmed in Fig.Q] where we compare the 
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FIG. 1: The CIMC ground-state energy estimates versus par- 
ticle number N for A/" m ax = 6 at unitarity using the renormal- 
ized Hamiltonian described in the text for different choices of 
the guiding wave function: PHFB (solid circles), PBCS (open 
squares), and HF (open diamonds). The statistical errors are 
smaller than the size of the symbols in all cases. 



CIMC ground-state energy estimates at the unitary limit 
for different choices of the guiding wave function: HF, 
PBCS and PHFB in a truncated model space of A/" max = 6 
oscillator shells. The energy estimates provided by the 
PHFB guiding wave function are much lower than those 
of the HF and PBCS wave functions for all values of the 
particle number N. 
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Ground state 


energy [huj] 


^Vmax 


CIMC 


Exact 


3 


6 


8.639(7) 


8.601 




7 


11.176(4) 


11.021 




8 


12.292(7) 


12.179 


8 


2 


2.000(6) 


2.000 




3 


4.391(4) 


4.279 




4 


5.208(9) 


5.138 



TABLE I: Comparison of the CIMC ground-state energies 
versus the exact ground-state energies for a CI model space 
of A/"max = 3 and TVmax = 8, and for several values of the 
particle number N. For the CIMC energies, the numbers in 
parenthesis denote the statistical error in the last significant 
digit. The jVmax = 3 energies were obtained using the oxbash 
code 26], while the TVmax = 8 energies were obtained using a 
serial version of the diagonalization code developed in Ref . \24. 

Next we investigate the quality of the PHFB wave 
function as a guiding wave function. To this end, we 
compare in Table U the CIMC energies using the PHFB 
as the guiding wave function with the exact results for 
smaller values of A/" m ax and N, for which exact numerical 
diagonalization is possible [2J, [26( . We find that for an 
even (odd) number of particles the CIMC energies are 
within 1% (2-3%) of the exact results. The accuracy of 
our method is comparable to that of the r-space fixed- 
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FIG. 2: The ground-state energy at unitarity for N < 20 par- 
ticles. The CIMC results for A/"max = 9 (solid circles) are com- 
pared with the results of Ref. _27J (open squares) and Ref. |2|| 
(open diamonds). The statistical errors are smaller than the 
size of the symbols. 

For larger systems exact numerical diagonalization is 
not possible. However r-space fixed-node GFMC re- 
sults are available [27], [28( . In Fig. [2] we compare the 
CIMC ground-state energies using A/ ma x = 9 with these 
fixed-node GFMC calculations for N < 20. The CIMC 
ground-state energies using A/" ma x = 9 oscillator shells 
are 0.5 — 3% higher than those obtained in Ref. (27|. The 
ground state energy estimates from Ref. [28| are typically 
slightly higher than those in Ref. [27] but with much larger 
statistical error. At A/" max = 9 the differences between the 
CIMC energies between successive values of Af max become 
smaller than their statistical errors. We have not carried 
out any A/" ma x — > oo extrapolations in this work. 

In Fig. [3] we compare the energy-staggering pairing 
gaps A calculated in the CIMC method with the pairing 
gaps obtained from the two r-space GFMC calculations 
mentioned above [27], [28| ■ We have also included the 
pairing gaps obtained from the density functional theory 
calculations of Ref. [29|. The energy-staggering pairing 
gap for an odd- A" system is defined by 



A = E (N)-~{E (N- 



1)+E {N + 1)} 



(17) 



where E (N) is the ground-state energy for N particles. 
Our results seem to be consistent with the results of 
Ref. i3 and Ref. [H, although the CIMC statistical error 
is much smaller than the statistical errors of the r-space 
GFMC calculations. The non-smooth behavior a,t N = 7 
is probably due to the shell closure at N = 8. 

For the results shown here we used ~ 2 x 10 5 indepen- 
dent configurations for the calculation of each energy. 
In most cases, the initial ensemble consisted of the non- 
interacting ground state. We also verified in a few test 
cases that the final result does not depend on the choice 
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FIG. 3: Energy-staggering pairing gaps versus number of 
atoms N. Our CIMC results for iV max = 9 (solid circles) 
are compared with similar gaps calculated in Ref. [27] (open 
squares), Ref. |28| (open diamonds) and Ref. 29j (open trian- 
gles). The vertical bars describe the statistical errors. The 
CIMC errors are smaller than the size of the symbols. 



of the initial ensemble. The calculations for W max < 4 
were carried out on a single processor and took up to 
a few hours. For 7V max > 4 the calculations were per- 
formed on a Linux cluster. The largest calculation with 
A/max = 9 and N — 20 took about 900 cpu hours. We 
choose 8t ~ 0.1 for these calculations. With this choice 
the number of time steps required to project the ground 
state was ~ 100 and the number of decorrelation time 
steps was < 5. We did not find any significant depen- 
dence of the optimal 5t, projection time or decorrelation 
time on Af s or N. 

The nominal scaling of the computational effort for 
the CIMC method with an AGP guiding wave function 
is - N 2 (N 8 - N) 2 xN 3 ~ N 5 (Af s - N) 2 where the factor 
N 2 (Af s — N) 2 comes from the maximal number of pos- 
sible non-zero matrix elements in Eq. ([5]) for a given n, 
while the factor N 3 is the computational effort required 
to calculate a pfaffian or a determinant of an N x N ma- 
trix. The pre-factor in this scaling can vary significantly 
depending on the type of interaction. For example, in 
the system described above with interactions between up 
and down spins only, there are only {N/2) 2 {(N S -N)/2} 2 
non-zero terms in Eq. © and we only have to calculate 
the determinant of an N/2 x N/2 matrix. Thus, the 
computational effort in this case is about two orders of 
magnitude smaller than in the most general case. 

It is interesting to compare our method with 
the constrained-path Monte Carlo [30] and the full 
configuration-interaction quantum Monte Carlo [3l| 
methods. Both of these methods perform a stochas- 



tic projection of the ground-state wave function in Fock 
space similar to our method. In the constrained-path 
Monte Carlo method, the Hamiltonian is 'linearized' with 
the help of the Hubbard-Stratanovich transformation and 
the random walk is carried out in the continuous space 
of the auxiliary fields. The sign problem is circumvented 
with the help of a guiding wave function. However, 
the mixed energy estimate in this method is not an up- 
per bound on the true ground-state energy. In the full 
configuration-interaction quantum Monte Carlo method, 
the random walk is carried out in the discrete configura- 
tion space as in our method. No guiding wave function 
is used in this method, and the sign problem is mitigated 
using a configuration-annihilation algorithm. However, 
the computational effort in this method scales with the 
size of the many-particle space (i.e., it is exponential in 
Af s /N), although it is a fraction of the computational 
effort involved in a direct diagonalization. 

In conclusion, we have introduced a novel CIMC algo- 
rithm that provides an upper bound for the ground-state 
energy of a finite fermionic system in the CI approach. 
We argue that the AGP class of wave functions provides a 
good choice for guiding wave functions. We demonstrate 
CIMC for the trapped cold atom Fermi gas condensate 
at unitarity using the PHFB wave function as a guid- 
ing wave function. For a small number of particles and 
sufficiently small number of single-particle orbitals, we 
find that the CIMC ground-state energies are within a 
few percent of the exact results. Our CIMC results for 
number of particles N < 20 are consistent with previous 
coordinate-space GFMC calculations. 

We emphasize that the CIMC method described here 
is quite general and can be applied to other fermionic 
systems such as cold atoms in a deformed trap and for 
finite nuclei. It can also be used to calculate properties 
other than the ground-state energy, such as the average 
occupation numbers. It will be interesting to determine 
whether the energy upper bounds can be improved by 
using other choices for the geminals in the AGP wave 
function, or by using other classes of wave functions al- 
together. 
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